Fluctuations and correlations in population models with age structure 
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We study the population profile in a simple discrete time model of population dynamics. Our model, 
which is closely related to certain "bit-string" models of evolution, incorporates competition for 
resources via a population dependent death probability, as well as a variable reproduction probability 
for each individual as a function of age. We first solve for the steady-state of the model in mean field 
theory, before developing analytic techniques to compute Gaussian fluctuation corrections around 
the mean field fixed point. Our computations are found to be in good agreement with Monte-Carlo 
simulations. Finally we discuss how similar methods may be applied to fluctuations in continuous 
time population models. 
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The problem of population dynamics has attracted 
enormous interest over many years (for some introduc- 
tions and recent applications see Refs. [|l]-f|). Beginning 
with simple logistic growth models [0 , a tremendous va- 
riety of systems have been studied displaying a diverse 
range of behavior, varying from stable fixed points to 
strange attractors. Broadly speaking these models split 
naturally into two categories: those using continuous and 
those using discrete time. The simplest discrete time 
models describe species where there is no overlap be- 
tween successive generations, leading to difference equa- 
tions of the form N(t + 1) = g[N(t)], where N(t) is 
the total population at time t. However these models 
may easily be generalized to species with multiple dis- 
crete age generations (for example: eggs, larvae, adults), 
where one or more generations may be present simul- 
taneously. Instead of a single variable N, information 
about the age distribution is now carried in a "vector" 
n(t) = {no(t),ni(t), . . . ,no(i)}, where n a (t) is the num- 
ber of individuals of age a at time t. Note that D is the 
maximum age, while no stands for the number of "new- 
borns" . Beginning with the pioneering work of Leslie @] 
models of this type have been extensively analyzed. The 
simplest Leslie model is linear in n, so that the evolution 
equation is just n(t + 1) = An(t). Here, A is the Leslie 
matrix 
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where the elements f a are the fecundities (number of off- 
spring produced) of individuals of age a, and the v a are 
Verhulst factors (the fraction of individuals of age a who 
survive to become age a + 1). Though the original Leslie 
model had all the {fa] and {v a ] as constants, general- 
izations have since been made to n-dependent factors 



0-12 1, so that the evolution dynamics becomes inher- 



ently non-linear. For example, n-dependent Verhulst 
factors are often used to mimic competition for finite re- 
sources. 

However one deficiency of the models discussed hith- 
erto is that they are deterministic. Real population sys- 
tems are of course affected by random fluctuations, com- 
ing from the environment and/or from the intrinsic dy- 
namics of the birth/death processes. Such stochastic 
Leslie models have also been investigated |||l3],[l4| , how- 
ever only for cases where the birth/death probabilities 
were independent of the population vector n. To date 
no information has been available regarding the more re- 
alistic case of fluctuations in stochastic age structured 
models with population dependent birth/death probabili- 
ties [^4). This is the situation we will study in this letter. 

Population models have also been intensively stud- 
ied by physicists in recent years in the context of so- 
called "bit-string" models of evolution M . These models 
are based on the mutation accumulation hypothesis [pi, 
which assumes that during the aging process each in- 
dividual accumulates exclusively late-acting deleterious 
genetic mutations. In "bit-string" models the genome 
of a particular species is encoded as a series of '0's and 
'l's (deleterious mutations), and as an individual ages 
the bits (genes) are activated one by one. When the 
accumulated sum of bad genes reaches a certain num- 
ber the individual dies (although death may also occur 
at younger ages due to Verhulst competition). Note that 
the "bit-strings" of the offspring may differ from those of 
the parent due to additional beneficial ('1'— >'0') or dele- 
terious ('0'— mutations. This type of model is clearly 
well suited to efficient computer simulation In its 

simplest case, where individuals die after the first delete- 
rious mutation, "bit-string" models simply correspond to 
multiple genome population models with age structure, 
where the different genomes can be distinguished by dif- 
ferent maximum ages. Deterministic versions of some of 
these models have already been treated analytically 0. 
However, an analysis of the important role played by fluc- 
tuations has so far been lacking. Our calculations form 
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the first step towards filling this gap. 

We begin our analysis by defining our discrete time 
population model. For simplicity, we consider only a 
single species reproducing asexually, without mutations. 
Thus, our system can be described by a single vector 
n. At each time step we compute the Verhulst factor 
V(n) and let each individual survive with probability V. 
After this "pruning" , each of the remaining individuals 
of age a may give birth to F a offspring with probability 
r a . At this point the remaining population is aged by 
one time step, with the exception of the new offspring 
who make up hq. Individuals who exceed the maximum 
age D then die immediately and are removed from the 
system. Since our model allows for a variable reproduc- 
tion probability function of age, features such as 
puberty and menopause can be naturally incorporated. 
However, we do assume that reproductive individuals of 
the same age produce an identical number (F a ) of off- 
springs. Note that we are not restricting ourselves to 
specific forms for V, r a , or F a beyond some general fea- 
tures, so that this analysis may easily be applied to the 
various "bit-string" models (|J^] . These general features 
include V,r a £ [0,1], since they represent probabilities. 
Furthermore, we assume that V depends only on the to- 
tal population N = ^2 a n a , via the ratio N/Nq, where 
No represents a characteristic population size that the 
resources can support. Also, to be reasonable, V is as- 
sumed to be monotonically decreasing with N . For com- 
parisons with simulations, we use V — sq exp(— N/Nq), 
where sq is another constant, a form frequently chosen in 
the biology literature [fij. In contrast, the algebraic form 
V = l — (N/No) is the favorite in the recent physics litera- 
ture |(||7j . We prefer the exponential form, since absolute 
cut-offs seem unrealistic in a real population system. 

As mentioned above, deterministic versions of this 
model have been studied |lo|-[l2]|; in particular the pop- 
ulation dynamics of a semivoltine species studied in 
Rcf. is quite similar to a D — 1 version of our model. 
Our goal is to go beyond these deterministic treatments 
and analyze the fluctuations and correlations in this sys- 
tem. Therefore, we need to consider P(n, t), the prob- 
ability of finding the population with a particular dis- 
tribution n at time t. Its evolution obeys the master 
equation 
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Note that the no+i is just a "temporary" variable, which 
keeps track of the number of mu's who survive the Ver- 
hulst "pruning" so that they can give birth before dying 



from old age. Let us emphasize that this equation is ac- 
tually quite complex, since V is a function of the total 
population. Multiplying (^) by n c and summing over all 
the other indices, we obtain 
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{no)t+i = F d r d (V{N) 
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where (») t denotes the average of • over P(n, t). These 
equations are exact. However, due to the presence of N 
through V, all moments of P may be coupled together. 
The mean field (MF) approximation consists of replacing 
the higher order moments by appropriate products of the 
first moment. Hence we find 



Wt+1 



\ n o) t +i 



\ MF 



(a > 0), (5) 



) t MF ) 



E P ; \ MF 

F d r d (n d ) t . 



(G) 



where, to be clear, we have written the explicit expres- 
sion for N . These non-linear equations are known to 
contain a rich variety of behavior, depending on the de- 
tails of F c ,r c , and V. For example, if ^ c F c r c < 1, the 
reproductive rates are too low and the population will 
eventually die out. On the other hand, if the reproduc- 
tive rates are large enough, the population will display 
period doubling bifurcations and chaos (l5). Let us fo- 
cus on the "moderate" range, so that a simple non-zero 
steady-state exists. In that case these equations are eas- 
ily solved to give 
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where z is the unique, positive, real root of the equation 
J2 C FcT c z c+1 = 1, and N(z) is the steady state total pop- 
ulation, given by the value that satisfies V(N) = z. Our 
objective is to go beyond such well-known mean field so- 
lutions and investigate fluctuations and correlations, i.e., 
the second moments of P. Thus, we multiply Eq. (|^) by 
n a ni, and sum over all the other indices. For a > 0, b > 0, 
we have 

(n a n b ) t+1 = (V 2 n a -inb-i) t + S ab (V(l - V)7i a -i) t , (8) 
(n a n ) t+1 = ^2 F c r c (V 2 n a ^ 1 n c ) t + 

C 

+F a -ir a - 1 {V(l-V)n a -i) t , (9) 
(non ) t+1 = £ F c r c F d r d (V 2 n c n d ) t + 
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+ J2[ F c( r c{Vn c ) t -r 2 c (V 2 n c ) t )]. (10) 

c 

Note that, like Eqs. ([|) and ([IJ), these are exact. Assum- 
ing No 3> 1, and that the system is well away from "criti- 
cal" points (e.g., bifurcations and the survival/extinction 
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transition), it is reasonable to postulate a Gaussian dis- 
tribution for P* (n) with width of 0(\fN~o). Rewriting 
Eqs. and (|^|l0|) for (•)*, we then have a closed set 

of equations for the first and second moments. This ap- 
proach should form the first step in a systematic expan- 
sion of all quantities in decreasing powers of No- Fur- 
thermore, in the same spirit, we will let n/No assume 
continuous values. As a check, we will compare the re- 
sults from this approach with those from a Monte-Carlo 
simulation, for a simple case. 
Proceeding, let us write 
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where we expect the unknown (to be determined) pa- 
rameters n and G to be of O(No) and O(l), respectively. 
Note that we will integrate n from — oo — * oo rather than 
from — > oo, a simplification which will introduce only 
negligible errors of 0(exp[— No])- Averages can now be 
computed using 
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Note that the second term in Eq. fll2j ) is expected to be 
suppressed by a factor of 0(1 /No) compared to the first 
term f(n a )- Hence the right hand side of Eq. Jl^ ) is actu- 
ally an expansion in powers oH/Nq. This ordering allows 
us to set up a systematic perturbation theory, which can 
be pushed to higher orders if desired. 

From now on we drop the bars for clarity (n a — > 
n a ). Defining £ = J2c n c/^o an( i V = dV/d£, we have 
OV/dn c = V'(0/N and d 2 V/dn c dn d = V"(0/N 2 , in- 
dependent of c or d. Applying Eq. ( |l^ ) to Eq. (||) gives 
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An equation for no can be similarly derived. With the 
assumptions n ~ O(Nq) and G ~ 0(1), the latter two 
terms in Eq. ( [l3"| ) represent 0(1 /No) corrections to the 
mean field results, while the remaining (lowest order) 
pieces make up the mean field equation (jq). Writing a 
perturbative expansion: n a = + n a + • ■ ■ (with 
n a k) ~ 0(N^ k )), we see that n a 0) is given by Eq. (@), 
while the first order result is 
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where I is the unit matrix and S is the stability matrix 
associated with the mean field (zeroth order) stationary 
solution. Note that both S and U need to be evaluated 
at zeroth order only. With our assumptions about F a , r a , 
and V(N), the eigenvalues of the stability matrix § usu- 
ally lie within the unit circle, implying that our mean 
field solution is stable. However, for sufficiently high 
reproductive rates, perturbations with 5n a oc n a (i.e., 
populations with the same relative age distribution, but 
with with different total numbers of individuals) can have 
eigenvalues of less than — 1. This is the signal of a pe- 
riod doubling bifurcation, leading to the breakdown of 
our Gaussian perturbation expansion (see also below). 

Applying the same analysis to the second moments, we 
obtain, after some lengthy algebra, 
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Again, all quantities need to be evaluated only at the 
zeroth order, so that, e.g., V is just z. In compact form 
this equation can be written as G — SG § T = K, which 
may be solved by series 

G = K + SKS T + § 2 K(§ T ) 2 + ... = ^§"K(§ T )". 

n 

(19) 

Since the eigenvalues and eigenvectors of S are known, let 
us write § = MEM -1 , where E is in Jordan form, with 
the eigenvalues on the diagonal, and M is the matrix 
(with its columns) composed of the corresponding right 
eigenvectors. Note that M is not necessarily orthogonal 
or unitary. If we define G = M _1 G (M T ) _1 and K = 
M _1 IK (M" 1 ) 7 , then it is straightforward to show that 
G = ^ n E"KE n . For simplicity, let us focus on the case 
where E is diagonal. Then the sum is easily performed, 
so that 



(20) 



Gab = K ab / (1 - e a e fc ) (no sum), 



where the {e} are the eigenvalues. Since G = I 
we can directly obtain the matrix G and with it all the 
information about the Gaussian probability distribution 
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(pi]). The explicit formula for computing G is our prin- 
cipal result. Given a particular form of V(N) and repro- 
ductive parameters r a ,F a , we can compute G and find 
the fluctuations in, as well as the correlations between, 
the populations of various ages. 

The result @ contains a further appealing feature: 
the signal of bifurcation. From stability analysis, 
we know that period doubling emerges when the eigen- 
value associated with Sn a oc n a reaches — 1. Examining 
Eq. (|2"c|), we see that it is precisely this feature which 
signals the breakdown of the Gaussian approximation. 
Furthermore, in many studies of, e.g., the Penna "bit- 
string" model |^| , menopause sets in before death, so 
that E is not diagonal. Then the final expression for 
G will be slightly more complicated, although the above 
conclusions will remain qualitatively unchanged. 

To check the above analysis, we study the simplest pos- 
sible case: a 2 age system (i.e. D = 1), with r a = F a = 1. 
Choosing the exponential form for V with sq = 1 and 
N = 100, mean field theory yields = 157.4 and 
rii = 119.3. Performing our analysis, we arrive at the 
first order corrections to no & n d ni, the fluctuations in 
the populations of each age, and the correlation between 
the populations of the two ages. The results are listed 
in Table EjL alongside those from Monte-Carlo simulations 
. The agreement is excellent, validating our approach. 

Note that the corrections n a /ria are less than 1%, vin- 
dicating our assertion that they should be 0(1/Nq). 

Up to this point we have been considering models with 
discrete time steps. However it is perfectly possible, and 
sometimes more appropriate biologically, to analyze mod- 
els in continuous time |^],^| . Let us conclude with some 
brief remarks about fluctuations in this context. A suit- 
able equation for the mean field population dynamics is 



TABLE I. Comparison of results for the 2 age model. 
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with boundary conditions for birth at x = and certain 
death at x = D. However the birth/death/aging pro- 
cesses giving rise to Eq. (^JJ) can also be written as a bal- 
listic reaction model on a discrete spatial lattice but with 
continuous time. As is well known fi"7j| , starting from the 
corresponding microscopic lattice master equation, tech- 
niques now exist to map this model onto a field theory 
in continuous space-time. The ensuing action can be 
recast as a Langevin equation, with the result being Eq. 
(pl|), but with extra multiplicative noise terms. The form 
of these noise terms would then be completely specified, 
without any ad-hoc guesses. Unfortunately the field- 
theoretic action is rather awkward, due to the presence of 
non-local interactions and non-local, multiplicative noise 
(from fluctuations in the birth process at x = 0). How- 
ever, simplifications occur if we are interested only in the 
simple, non-zero steady state, where expansions about 
the mean field solution should be adequate. In this case 
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the leading noise terms enter additively, so that a pertur- 
bation theory analogous to the above approach can be 
set up. 

In conclusion, we have developed analytic techniques 
for dealing with fluctuation effects in a general class of 
population models with age structure. The results we 
have presented also form a first step towards an improved 
analytic understanding of the "bit-string" models of evo- 
lution. Finally, it would be interesting to perform further 
investigations near the bifurcation point, since interest- 
ing collective behavior can be expected in that region. 
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